{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "dee6bee1",
   "metadata": {},
   "source": [
    "# successful try od fiducial features"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea30f0d3",
   "metadata": {},
   "source": [
    "## tompkins do automatic segmentation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28d69f0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.signal import butter, filtfilt, find_peaks\n",
    "\n",
    "def butter_bandpass(lowcut, highcut, fs, order=5):\n",
    "    nyq = 0.5 * fs\n",
    "    low = lowcut / nyq\n",
    "    high = highcut / nyq\n",
    "    b, a = butter(order, [low, high], btype='band')\n",
    "    return b, a\n",
    "\n",
    "def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):\n",
    "    b, a = butter_bandpass(lowcut, highcut, fs, order=order)\n",
    "    y = filtfilt(b, a, data)\n",
    "    return y\n",
    "\n",
    "def pan_tompkins(ecg_signal, fs=1000):\n",
    "    # Bandpass filter the ECG signal\n",
    "    ecg_signal_filtered = butter_bandpass_filter(ecg_signal, 5, 15, fs)\n",
    "\n",
    "    # Differentiate the signal\n",
    "    diff_signal = np.diff(ecg_signal_filtered)\n",
    "\n",
    "    # Square the signal\n",
    "    squared_signal = diff_signal ** 2\n",
    "\n",
    "    # Integrate the signal\n",
    "    integrated_signal = np.convolve(squared_signal, np.ones(fs // 8))[:len(ecg_signal)]\n",
    "# ######\n",
    "    # Find the QRS complex peaks using a threshold-based approach\n",
    "    qrs_peaks, _ = find_peaks(integrated_signal, height=0.35*np.max(integrated_signal), distance=0.2*fs)\n",
    "\n",
    "    # Calculate the fiducial features\n",
    "    qrs_onsets = []\n",
    "    qrs_offsets = []\n",
    "    p_peaks = []\n",
    "    t_peaks = []\n",
    "    rr_interval_low = int(fs * 0.2)\n",
    "    rr_interval_high = int(fs * 1)\n",
    "    \n",
    "    for qrs_peak in qrs_peaks:\n",
    "        # Find the QRS onset\n",
    "        # Find the QRS onset and offset\n",
    "        qrs_onset = qrs_peak - int(rr_interval_low * 0.5)\n",
    "        qrs_offset = qrs_peak + int(rr_interval_high * 0.5)\n",
    "        qrs_onsets.append(qrs_onset)\n",
    "        qrs_offsets.append(qrs_offset)\n",
    "\n",
    "#         # Find the QRS onset\n",
    "#         for i in range(qrs_peak, 0, -1):\n",
    "#             if integrated_signal[i] < integrated_signal[qrs_peak] * 0.15:\n",
    "#                 qrs_onsets.append(i)\n",
    "#                 break\n",
    "#         # Find the QRS offset\n",
    "#         for i in range(qrs_peak, len(ecg_signal)):\n",
    "#             if integrated_signal[i] < integrated_signal[qrs_peak] * 0.2:\n",
    "#                 qrs_offsets.append(i)\n",
    "#                 break\n",
    "                \n",
    "                \n",
    "        # Find the P wave peak\n",
    "        p_peak, _ = find_peaks(ecg_signal_filtered[qrs_onsets[-1]:qrs_peak], height=0.1*np.max(ecg_signal_filtered), distance=0.1*fs)\n",
    "        if len(p_peak) > 0:\n",
    "            p_peaks.append(p_peak[-1] + qrs_onsets[-1])\n",
    "        else:\n",
    "            p_peaks.append(None)\n",
    "        # Find the T wave peak\n",
    "        t_peak, _ = find_peaks(-ecg_signal_filtered[qrs_peak:qrs_offsets[-1]], height=-0.1*np.max(ecg_signal_filtered), distance=0.2*fs)\n",
    "        if len(t_peak) > 0:\n",
    "            t_peaks.append(t_peak[0] + qrs_peak)\n",
    "        else:\n",
    "            t_peaks.append(None)\n",
    "\n",
    "    return qrs_peaks, qrs_onsets, qrs_offsets, p_peaks,t_peaks,integrated_signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "654a291e",
   "metadata": {},
   "outputs": [],
   "source": [
    "qrs_peaks, qrs_onsets, qrs_offsets, p_peaks,t_peaks,integrated_signal =pan_tompkins(sig3, fs=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5685464d",
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_array = []\n",
    "ecg_signal_filtered=filtered_signal.copy()\n",
    "for i in range(len(qrs_peaks)):\n",
    "    feature_vector = []\n",
    "    # Add QRS complex features\n",
    "    feature_vector.append(ecg_signal_filtered[qrs_peaks[i]])\n",
    "    feature_vector.append(integrated_signal[qrs_peaks[i]])\n",
    "    feature_vector.append(qrs_offsets[i] - qrs_onsets[i])\n",
    "    feature_vector.append(qrs_onsets[i] - qrs_onsets[0])\n",
    "    \n",
    "    if p_peaks[i] is not None:\n",
    "        feature_vector.append(ecg_signal_filtered[p_peaks[i]])\n",
    "        feature_vector.append(qrs_peaks[i] - p_peaks[i])\n",
    "    else:\n",
    "        feature_vector += [None, None]\n",
    "    # Add T wave features\n",
    "    if t_peaks[i] is not None:\n",
    "        feature_vector.append(ecg_signal_filtered[t_peaks[i]])\n",
    "        feature_vector.append(t_peaks[i] - qrs_peaks[i])\n",
    "    else:\n",
    "        feature_vector += [None, None]\n",
    "    feature_array.append(feature_vector)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f7dfbee",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(feature_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "826d25e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(feature_array[3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a554e3d5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "90e5df84",
   "metadata": {},
   "source": [
    "# another try to calculate fiducial features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2dca0ea7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.signal import butter, filtfilt\n",
    "\n",
    "def butter_bandpass(lowcut, highcut, fs, order=5):\n",
    "    nyq = 0.5 * fs\n",
    "    low = lowcut / nyq\n",
    "    high = highcut / nyq\n",
    "    b, a = butter(order, [low, high], btype='band')\n",
    "    return b, a\n",
    "\n",
    "def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):\n",
    "    b, a = butter_bandpass(lowcut, highcut, fs, order=order)\n",
    "    y = filtfilt(b, a, data)\n",
    "    return y\n",
    "\n",
    "def pan_tompkins(ecg_signal, fs=1000):\n",
    "    # Bandpass filter the ECG signal\n",
    "    ecg_signal_filtered = butter_bandpass_filter(ecg_signal, 2, 50, fs)\n",
    "\n",
    "    # Differentiate the signal\n",
    "    diff_signal = np.diff(ecg_signal_filtered)\n",
    "\n",
    "    # Square the signal\n",
    "    squared_signal = diff_signal ** 2\n",
    "\n",
    "    # Integrate the signal\n",
    "    integrated_signal = np.convolve(squared_signal, np.ones(fs // 8))[:len(ecg_signal)]\n",
    "\n",
    "    # Find the QRS complex peaks using a threshold-based approach\n",
    "    qrs_peaks, _ = find_peaks(integrated_signal, height=0.35*np.max(integrated_signal), distance=0.2*fs)\n",
    "\n",
    "    # Calculate the fiducial features\n",
    "    qrs_onsets = []\n",
    "    qrs_offsets = []\n",
    "    for qrs_peak in qrs_peaks:\n",
    "        # Find the QRS onset\n",
    "        for i in range(qrs_peak, 0, -1):\n",
    "            if integrated_signal[i] < integrated_signal[qrs_peak] * 0.15:\n",
    "                qrs_onsets.append(i)\n",
    "                break\n",
    "        # Find the QRS offset\n",
    "        for i in range(qrs_peak, len(ecg_signal)):\n",
    "            if integrated_signal[i] < integrated_signal[qrs_peak] * 0.2:\n",
    "                qrs_offsets.append(i)\n",
    "                break\n",
    "\n",
    "    return qrs_peaks, qrs_onsets, qrs_offsets \n",
    "\n",
    "def extract_fiducial_features(ecg_signal, fs=1000):\n",
    "    # Detect QRS complexes using the Pan-Tompkins algorithm\n",
    "    qrs_peaks, qrs_onsets, qrs_offsets = pan_tompkins(ecg_signal, fs=fs)\n",
    "\n",
    "    # Calculate other fiducial features, such as P and T wave peaks, if needed\n",
    "\n",
    "    return qrs_peaks, qrs_onsets, qrs_offsets\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1eb72cbe",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.signal as signal\n",
    "\n",
    "def detect_t_wave_peaks(ecg_signal, qrs_onsets, qrs_offsets, fs=1000):\n",
    "    t_wave_peaks = []\n",
    "    for i in range(len(qrs_onsets)):\n",
    "        qrs_complex = ecg_signal[qrs_onsets[i]:qrs_offsets[i]]\n",
    "        st_segment = ecg_signal[qrs_offsets[i]:qrs_offsets[i]+int(0.2*fs)]\n",
    "        st_mean = np.mean(st_segment)\n",
    "        st_std = np.std(st_segment)\n",
    "        threshold = st_mean + 0.5 * st_std  # adjust threshold as needed\n",
    "        t_wave_peak, _ = signal.find_peaks(ecg_signal[qrs_offsets[i]:qrs_offsets[i]+int(0.6*fs)], height=threshold)\n",
    "        if len(t_wave_peak) > 0:\n",
    "            t_wave_peak = t_wave_peak[0] + qrs_offsets[i]\n",
    "            t_wave_peaks.append(t_wave_peak)\n",
    "    return t_wave_peaks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a93c732b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def extract_features(ecg_signal, qrs_peaks, qrs_onsets, qrs_offsets, fs=1000):\n",
    "    # Calculate heart rate variability (HRV)\n",
    "    rr_intervals = np.diff(qrs_peaks) / fs\n",
    "    mean_rr = np.mean(rr_intervals)\n",
    "    std_rr = np.std(rr_intervals)\n",
    "\n",
    "    # Calculate QRS duration\n",
    "    qrs_durations = (qrs_offsets - qrs_onsets) / fs\n",
    "    mean_qrs_duration = np.mean(qrs_durations)\n",
    "    std_qrs_duration = np.std(qrs_durations)\n",
    "\n",
    "    # Calculate QT interval\n",
    "    t_wave_peaks = detect_t_wave_peaks(ecg_signal, qrs_onsets, qrs_offsets, fs)\n",
    "    qt_intervals = (t_wave_peaks - qrs_onsets) / fs\n",
    "    mean_qt_interval = np.mean(qt_intervals)\n",
    "    std_qt_interval = np.std(qt_intervals)\n",
    "\n",
    "    # Calculate ST segment deviation\n",
    "    st_segments = []\n",
    "    for i in range(len(qrs_peaks)-1):\n",
    "        start = int((qrs_peaks[i] + qrs_offsets[i])/2)  # mid-point of the QRS complex\n",
    "        end = start + int(0.2 * fs)  # 0.2 seconds after the mid-point\n",
    "        st_segments.append( ecg_signal[int(start):int(end)] - np.mean(ecg_signal[int(start-0.1*fs):int(start)]))\n",
    "  # deviation from baseline\n",
    "    st_segments_array = np.vstack(st_segments)\n",
    "    mean_st_deviation = np.mean(st_segments_array)\n",
    "    std_st_deviation = np.std(st_segments_array)\n",
    "\n",
    "    return {\n",
    "        'mean_rr': mean_rr,\n",
    "        'std_rr': std_rr,\n",
    "        'mean_qrs_duration': mean_qrs_duration,\n",
    "        'std_qrs_duration': std_qrs_duration,\n",
    "        'mean_qt_interval': mean_qt_interval,\n",
    "        'std_qt_interval': std_qt_interval,\n",
    "        'mean_st_deviation': mean_st_deviation,\n",
    "        'std_st_deviation': std_st_deviation,\n",
    "    }\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f5d79fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "qrs_peaks, qrs_onsets, qrs_offsets=extract_fiducial_features(filtered_signal, fs=1000)\n",
    "\n",
    "res=extract_features(filtered_signal,np.array(qrs_peaks), np.array(qrs_onsets), np.array(qrs_offsets), fs=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0c4a088",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "nlp",
   "language": "python",
   "name": "nlp"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
